!****************************************************************************
!
!  Code for Bianchi, Ottonello and Presno (2021) - Fiscal Stimulus under Sovereign Risk
! 
!  Updated: 12/05/2021
!
!****************************************************************************

program MAIN
    
USE GLOBAL
!USE FUNCALLER
USE OMP_LIB

implicit none

REAL(8)::uncstda
REAL(8)::dif,difaut,difq,time
INTEGER(4)::i,j,ii,jj
INTEGER(4)::i_b, i_a,ia0,ib0,iter,iterv
INTEGER::hz,ier,i_a0,i_b0(1)
REAL(8)::ones(na,nb),dos
DOUBLE PRECISION :: a_valor, b_valor,f_valor,indicator_external
REAL(8)::fvalue,auxval,func_gN,func,func_COMMIT2,func_COMMIT3, func_COMMIT3_LOWERBOUND,func_COMMIT3_UPPERBOUND,func_COMMIT3_SIZECUT,func_COMMIT3_debt,func_COMMIT3_yT
EXTERNAL  func_gN,func_COMMIT2 
REAL(8)::a_initial,b_initial,ftol
INTEGER::i_default_global
REAL(8)::beta0,chi0,icost0,psi10,psi20,wbar0,xval(6),xi(6,6)
REAL(8)::p0,cnval,gnval,fvalor

INTEGER::end_time,cr     

iterfun=0
choice_glob=1

!baseline:
bmin = -0.41d0          
bmax =  0.35d0

version_code=312
      
!        beta          chi       icost        psi1         psi2        wbar  !
xval = [0.907d0,      0.02d0,    0.0d0,       0.3772d0,    2.42d0,     3.068d0 ]

call system_clock(count_rate=cr)
rate_time = REAL(cr)
call system_clock(count=start_time)

call quadrature

uncstda = sigmaa/DSQRT(1d0-rhoa**2d0)    
amin = mua/(1d0-rhoa) - ndevs*uncstda
amax = mua/(1d0-rhoa) + ndevs*uncstda

DO i=1,ns
    agrid(i)= amin + ((REAL(i)-1d0)/(REAL(ns)-1d0))*(amax - amin)
ENDDO    

OPEN(1,FILE='output//agrid.txt')
OPEN(2,FILE='graphs//agrid.txt')

DO i = 1,na
    WRITE(1,*) agrid(i)
    WRITE(2,*) agrid(i)
ENDDO
    
CLOSE(1)
CLOSE(2)

DO i=1,nn
    gngrid(i) = gnmin + ((REAL(i)-1d0)/(REAL(nn)-1d0))*(gnmax - gnmin)
ENDDO 

DO i=1,nn
    hgrid(i) = hmin + ((REAL(i)-1d0)/(REAL(nn)-1d0))*(hmax - hmin)
ENDDO 

do i=1,ncdf
   cdfgrid(i) = cdfmin + (cdfmax - cdfmin) * (i-1) / (ncdf-1)
end do

dos = 2d+00
epsmin = -4d0*sigmaa
epsmax = 4d0*sigmaa

!grid when using Gauss-Hermite quaradture points 
do i=1,neps
   epsgrid(i) = SQRT(dos)*sigmaa* quad_x_hermite(i)
end do

kw   = omegac/(1d0-omegac)
payoff = 1d0
IF (choice_coup.EQ.1) payoff = lambda+(1d0-lambda)*coup

defaultgrid(1)=0d0
defaultgrid(2)=1d0

ftol=1E-3

xi(:,:)=0d0
FORALL (i=1:6, j=1:6, j.EQ.i) xi(i,j)=1d0

fvalor = func(xval)
fvalor = func_COMMIT3(xval)
STOP

!compute welfare gains as function of cut size 
!fvalor= func_COMMIT3_SIZECUT(xval)

!compute welfare gains of blunt cuts as function of current yT
!fvalor=func_COMMIT3_YT(xval)
!compute welfare gains of blunt cuts as function of current debt
!fvalor=func_COMMIT3_debt(xval)
!
!To compute current adjusted gN that makes govt indifferent between repayment and default:
!(run func first to upload parameter values)
!as function of current yT
!fvalor=func(xval)
!CALL WELFARE_GAIN_DEF_YT

!as function of current debt
!fvalor=func_gN(xval)
!CALL WELFARE_GAIN_DEF_DEBT
!STOP
!
!fvalor=func_gN(xval)
!!fvalor=func_COMMIT3(xval)
!compute welfare gains state-contingent cuts as function of upper bound for total income 
!fvalor=func_COMMIT3_UPPERBOUND(xval)
!compute welfare gains state-contingent cuts as function of lower bound for total income 
!fvalor=func_COMMIT3_LOWERBOUND(xval)

call system_clock(count=end_time)      
PRINT*,'Time Elapsed =',(end_time-start_time)/rate_time

end program MAIN

